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Abstract. The nonlinear saturation of the tearing mode is revisited in slab 
geometry by taking into account higher-order harmonics in the outer solution. 
The general formalism for tackling this problem in the case of a vanishing current 
gradient at the resonant surface is derived. It is shown that, although the higher- 
order harmonics lead to corrections in the final saturation equation, they are of 
higher order in the perturbation parameter, which provides a formal proof that 
the standard one-harmonic approach is asymptotically correct. 
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1. Introduction 

The tearing mode [1] is a resistive magnetohydrodynamic (MHD) instability that 
resonates on magnetic surfaces where the mode's wave vector is perpendicular to the 
magnetic field. It leads to the formation of so-called magnetic islands around the 
resonant surface, resulting in a local change of magnetic topology and a greater radial 
transport on a length scale of the order of the islands' full width, w. Therefore, these 
structures are of great interest for nuclear fusion, and it is especially important to 
have a keen insight into their stability properties as well as their maximum achievable 
amplitude. 

While the linear theory of the tearing mode is rather well understood, the 
nonlinear one, pioneered by Rutherford thirty-five years ago [2], is taking longer to 
unravel. Recently, a series of works have derived solutions for the saturation of the 
tearing mode in the framework of simple physical models [3-8] , somewhat rekindling 
interest in this subject. Although rigorous in their mathematical details, they share a 
common feature based on an assumption first made by Rutherford, i.e. they neglect 
higher-order (poloidal) harmonics in the magnetic perturbation. However, this is not 
really justified a priori, since nonlinearities naturally couple all harmonics, and one 
may therefore question the validity of the afore-mentioned references. 

In this work, we investigate this problem using the simplest of physical models, 
namely that of reduced MHD in slab geometry in the so-called symmetric case (i.e. no 
current gradient at the resonant surface), which is the same as that used in [3,4]. We 
first derive the general solution to this problem and then explicitly solve the case with 
two harmonics for two different types of equilibrium. Last, we compare our results 
with those obtained from previous theories and draw conclusions. 
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2. Model equations 

We first introduce the following normalizations: 

t = r v t;x = Lx;y = Ly;J= J N J ; ip = p, J L 2 ip ; f = — <p, (1) 

Mo 

where t is the time variable, t v — p§L 2 jr\ is the resistive diffusion time, po is the 
permeability of free space, rj is the (uniform) resistivity, x and y are the radial and 
poloidal variables respectively, L is a characteristic radial length, Jjy = J eq {0), J 
(resp. J eq ) is the (resp. equilibrium) current density, ip is the magnetic flux function 
(i.e. B = B z e z + V x (ipe z ), where B is the magnetic field and e z is a unit vector 
perpendicular to the xy plane) and ip is the electric potential and plays the role of 
the (ion) stream function (i.e. v = e z x V<p, where v is the velocity field). Note that 
these normalizations are such that it is the equilibrium current density, and not the 
equilibrium magnetic field, that is normalized to unity at the resonant surface. We 
then use the reduced MHD equations in slab geometry [9] which, taking into account 
the normalizations while omitting the " ~ " for clarity, read: 

d t A ± ip + [<p, A ± <p] = S 2 [J,^] + ^Al^ (2) 

He 

d t l/} + [<p, 1p] = J eq - J (3) 

J = -A X V, (4) 

where 5* = vaLpo/v is the Lundquist number, Re = vaL/v is the Reynolds 
number, va = JnL\/ Ho/ P is the Alfven speed, v is the viscosity and p is the mass 
density (v and p are assumed to be constant). The Poisson brackets are given by 
[/, g] = d x fd y g — d x gd y f. Finally, the resonant surface is conveniently set at the 
origin by letting ip' eq (0) = 0, where the equilibrium magnetic flux function satisfies 

1p e q = — Jeq- 

3. Perturbed equilibrium 

Provided there are no equilibrium flows in the xy plane and given the fact that 
S,Re 3> 1, as is the case in present-day tokamak plasmas, the equation of motion 
([2]) simply yields: 

[J,tP]=0. (5) 
We then look for a perturbed equibrium of the form: 

V> ~ ipeq(x) + tl(t)Mx) cos [n X + &(«)], (6) 

n>l 

where we have defined x = ky, k being the mode's wave-number and, throughout 
this paper, '~' has the meaning 'equals plus higher order terms'. Note that, without 
loss of generality, it is possible to choose j3\ = 0, which will henceforth be the case. 
Substituting this expression into ([5]), the ip n satisfy: 



°eq 



i- C O 



/>Pn = 0- (7) 

L Yeq 

At this point, we should clarify the case of the n — component. Indeed, to be 
fully general, one may be tempted to include a term of the form 6Q(t)ipo(x) in ([6]) but, 
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because of its Poisson bracket nature, is automatically satisfied in order Sq and 
ipo( x ) remains unconstrained. Therefore, the n = component can only originate in 
Ohm's law's quasilinear terms, and one can show that this leads to an order ^n>i 
change to the equilibrium, which, as was already noted in [2] , can be neglected in (J5J . It 
has also been shown that the magnetic island itself may result in an n = component, 
but the latter has no impact whatsoever on the saturated island width [6,7]. Since 
the main interest of the present paper is the saturation of the tearing mode, we shall 
ignore the n — component altogether for the sake of clarity. 

In the following, we suppose that the first harmonic dominates the others at all 
times, i.e. that 8\ > 8 n> \. This assumption certainly makes sense in the linear regime, 
since the first harmonic is the most unstable one [1], but it also has to be consistent 
with the final saturation result, which will be shown to be the case later on. Based 
on the work already done in e.g. [6,7], one can then infer from ([3]) that there is a 
boundary layer of width 5\ centered on x = and that one thus has to resort to the 
technique of asymptotic matching. Writing J eq ~ 1 + a 2 x 2 , choosing ip eq {0) = for 
convenience and solving J7]| using Froebenius' method [10], we derive the following 
expansion for the outer solution: 

i! 

2 



Co 



E 

n>l 



a n cos (nx + Pn) 



A' 



S-i |f | 2^ a n ~y cos (nx + (3 n 

n>\ 



(8) 



51 



0,2 



n>l 



«2 



(nkf 



f 2 cos (nx + n ) 



where we have defined (out = — ot-n = (o~n/Si) 2 and the inner variable f = x/Si 
in order to make the ordering in the small parameter Si explicit. As to the linear 
stability parameters, A^, they are related to the logarithmic jump of the ip n around 
the resonant surface [1], namely 



a: 



lim -W n (e) 

e— »0+ 



<(-£)]/^n(0), 



(9) 



and depend on both the equilibrium current density profile and the boundary 
conditions. 

The solution ([5]) holds in the so-called outer region of the plasma where ideal 
MHD is a good approximation to our set of equations, but breaks down around the 
resonant surface where resistivity has to be taken into account. The solution that is 
valid in this resistive boundary layer is the inner one which we derive in the next 
section. It utlimately has to be matched to © that, in effect, is analogous to a 
boundary condition at infinity (i.e. |f| — + oo). At this stage, it is perhaps worthwhile 
to stress once more that it is the n > 1 terms in the sums appearing in §8§ that were 
neglected in previous works. 
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4. Solution in the inner region 

4--1- Inner equations 

The inner equations are obtained by re- writing with respect to the inner 

variable £: 

[C, J] = o (to) 
- 6fd t c + SxdM^c - 20 + fc<5i[C, <p] ~ i + a 2 sje - J (n) 
j = a|c + * a «?3fr (12) 

where the Poisson brackets are now taken with respect to the (£, %) variables. Note 
that IjlOp is valid in the nonlinear regime only, where the island is supposed to be 
larger than the resistive and visco-resistive layer widths [2]. This set of equations 
is then classically solved using perturbation expansions in powers of Si, i.e. writing 
C = E(>o S 'iCi, V = E;>o ^i- and J = E;>o s i J i- 



4.2. Order 

Ohm's law simply gives Jo = 1, and the integration of (fT2|) together with the matching 
condition provided by ([5]) at this order yields: 

Co =?/2-J2<x n cos (n X + /3n)- (13) 

n>l 



4.3. Order 5 1 

Equation (|10p implies that the order 1 component of J is a function of Co only on 
either side of the resonant surface, i.e. J\(£,x) = ii(Coji)j where we have defined 
± = sign(£). It is therefore easier to work in (Co,X!±) variables, which will be the 
case in the following. Ohm's law then yields: 

2 ^2 dtS n y/a^cos (nx + fin) + k£ <9 x <^ | Co = -jl(Co5 ±)- (14) 

In order to solve this equation, it is convenient to define, for any function /, its flux 
average (/) as: 

f dxf{z,x\±)/£ if2>C sep = Co(0,7r) 
< J ~* r xo , (15) 

1/2 £ o- / d X f(z, X ;<r)/t if^<Cse P 

a=± Xo 

where xo € [0, if] is the turning point of the corresponding flux surface, i.e. it satisfies 
Co(0,xo) = z : an d, in this expression, £ has to be taken as a function of (z, x;±), i.e. 
C = ±v2[z + En>i a ™ cos i n X + Pn)] 1 / 2 - With this definition in mind, it is then easy 
to show that the solution to ([14]) reads: 

h =-2^d t <W^(cos(nx + /9„))/<l>. (16) 

n>l 
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44. Order 5\ 

The calculation is similar to that of the previous section except that we now neglect 
all terms dependant on <9t<5i since they would only lead to higher order corrections in 
the final result. Consequently, ([3]) gives: 

- ^2 a n d t (3 n sin (nx + n ) + fc £ <9 x <^ol Co = a 2 £ - j 2 (Co; ±), (17) 
n>l 

which, after applying the bracket (■ • •) operator on both sides, yields: 



32 



«2<C 2 ) + X! a ndt/3 n {sin {nx + /?„)) 



/<!>• (18) 



We stop the inner calculation here since it is the lowest relevant order. Indeed, 
we shall see in the next section that the asymptotic matching conditions on the A' n 
terms in (jHJ already provide a saturation theory for the tearing mode at that order. 

5. Asymptotic matching conditions 

When taking the asymptotic expansion of the inner solution derived in Section [¥] for 
|£| — > 00, one shows that the matching with the outer solution given by © provides 
the following conditions (see, e.g., [7] for more details on the asymptotic matching 
procedure): 

~ j£(cos(m X + [3 m ))E(( )~a m A' m (19) 

K JUin (1) 
+ 00 i> 

7 ^(sin(mx + A n ))£(Co)~0, (20) 

where ( min = ( (0, 0) and E(( ) is given by: 

E = ji + Si ^2 a n (2a 2 cos (nx + (3„) + d t (3 n sin (nx + (3 n )) ■ (21) 
n>l 

Since we want to focus on the saturation of the mode, it is possible to simplify this set 
of equations. Indeed, (|20| implies that, at saturation, f3 n € {0, 7r} for all n > 1 (recall 
that (3\ = by assumption). Therefore, we can do away with the /3 n 's altogether 
provided we allow the a„'s to be either positive or negative (except for ct\ that is 
always equal to 1). If we do this, ([2H)l is automatically met and, letting dt — 0, (fH?)) 
can be re-written as: 

a m A' m + a 2<?l ^ Qnfmn({aj}) ~ 0, (22) 
n>l 

where we have defined 

F mn ({a.i\) = -\ — r (cosmx) (cosnx) • (23) 

Note that, trivially, F mn = F„ m . The set of equations given by (|22p provides a general 
theory for the nonlinear saturation of the full island width w s = 4c?! (2Zn>o a 2n+i) 1 ^ 2 
when taking into account any number of harmonics in the outer solution (the island 
width is defined as the width of the separatrix at the O-point, where the separatrix is 
given by the equation Co(£,x) = Csep = Co(0,tt) = E™>i(-l) n+1 «n)- 
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Although it is not possible to solve it analytically, we make an important comment 
about this result. Indeed, if we relax Rutherford's original assumption (namely 
— A' n > 1 for i! > 1) [2], and instead (more reasonably) assume A' n>1 to be of order 
one, then (|2"2")1 implies a n >i = 0(a 2 5\/ A' n ) = 0(5i). Consequently, the higher-order 
harmonics only lead to a higher order correction in the saturation equation, i.e.: 

A; + a 2 5 1 { F n + 5^ a n F ln J = Ai + oa^-Fii + 0{5\) ~ 0, (24) 

V n>l J 

showing that the standard one-harmonic approach is asymptotically correct in the 
limit of vanishing Si. Incidentally, the fact that a n >\ — 0(8i) also ensures that our 
original assumption, namely that the first harmonic dominates over the higher ones, 
is consistent. In order to demonstrate this result quantitatively, we now illustrate our 
theory through the simplest case of two harmonics and compare the outcome with 
that of [3,4]. 



6. Corrections due to the second harmonic for two specific equilibria 



Taking into account the first two harmonics only, (|22[) gives: 

4A'j +a 2 w s (F 11 + aF 12 ) =0 (25) 
4aA' 2 +a 2 w s (F 2 i + aF 22 ) = 0, (26) 

where now, since we truncate after the second harmonic, w s = A8\ and we have 
written a = a 2 for the sake of clarity. Since a is expected to be small, it makes sense 
to Taylor expand the F coefficients to first order as F mn ~ A mn + aB mn . It is then 
straightforward to solve for w s and a: 

-4A; _ A 12 A[ 

Ws ~ a 2 [A n + Q (B U + A 12 )} 5 a " A' 2 A n - A[(B 12 + A 22 ) ' (27) 

where A n ~ 3.29, B n ~ 2.22, A 12 ~ 0.34, B 12 ~ 2.03 and A 22 ~ 1.59 have been 
computed numerically. 

These expressions can be evaluated provided A' 1; A' 2 and a 2 are given. To this 
end, we consider two specific equilibria often used in the literature: the cosh (i.e. 
ipeq — 1/(2 cosh 2 x) [11], where the 1/2 factor has been included to comply with our 
normalization scheme) and sheet pinch (i.e. ip e q = — log (cosh a:) [12]) equilibria. The 
first one has a 2 = —4 and the other a 2 = —1, while the A^ are given by: 



2(5-n 2 fc 2 )(3 + n 2 fc 2 ) 
n 2 k 2 V 4 + n 2 k 2 



and Ai 



nk 



nk 



(28) 




Aw(%) 



Figure 1. Plot of Aw ( ) and a ( ) as a function of for the cosh (a) 

and sheet pinch (ft) equilibria. 
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respectively. The results so-obtained are then compared with the one given by 
standard theory, namely w s = — 4A' 1 /a 2 A 11 [3,4], and the outcome is shown in figure[TJ 
where we have plotted Aw = (w s — w s )/w s and a for both equilibria. We see that 
the island is found to be slightly bigger than was predicted by previous theory, but 
the corrections are very small, namely up to 5% (resp. 3%) larger for the cosh (resp. 
sheet pinch) equilibrium with < 5 (resp. < 1). Of course, the greater A^, 
the bigger these corrections can get, and, indeed, they reach up to 10% (resp. 15%). 
However the constant- ?/> approximation breaks down for too large a A[, so that the 
ranges shown in figure [1] are, in effect, appropriate as far as our asymptotic matching 
theory is concerned. Therefore, we conclude that the standard theory of [3,4] gives a 
correct result in its region of validity, which was expected since it has been confirmed 
numerically in [13]. 

We now make one last remark concerning these results. One may naturally wonder 
whether the neglect of the third and higher harmonics is justified, since (f!Z!?)) only 
implies that all n > 1 harmonics are 0{8\/ A' n ). To test that this is the case, we have 
run a full numerical simulation for the cosh equilibrium with the pseudospectral code 
already used in [13], taking A^ = 1.0. The results are shown in figureEJ where we see 
that the second harmonic always has a greater amplitude than the third and fourth 
(as well as the higher-order ones, not shown). The reason for that is twofold: first, 
the nonlinear coupling to the first harmonic can be inferred to grow weaker for higher- 
order ones, and second, since A' n oc n for large n, the amplitude of the higher-order 
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Figure 2. (colour online). Plot of the first four harmonics' amplitude versus time 
for the cosh equilibrium with A' 1 =1.0. 
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harmonics scales (at most) as 1/n. Incidentally, we also see in figure [2] that the first 
harmonic largely dominates at all times, which gives further evidence that our original 
assumption (namely Si > S n> i) is correct, at least in the regime (A^ not too large) 
we investigate. 

7. Conclusion 

In this paper, we have examined the impact of higher-order harmonics on the 
nonlinear evolution of the tearing mode in (symmetric) slab geometry. We have 
provided a general set of dynamical equations as well as a compact formula for 
the saturation of the mode. We have shown that the contribution due to the 
higher-order harmonics in the saturation equation led to a higher order correction 
in the perturbation parameter. Hence, we have justified the standard one-harmonic 
calculation as asymptotically valid. To make this claim more tangible, we have 
computed numerically the contribution due to the second harmonic for two types 
of equilibrium and shown that it indeed led to small corrections in a relevant range of 
values for A[. Therefore, we believe this work usefully complements previous results 
in the theory of tearing mode saturation. 
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